Optimization by thermal cycling 
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An optimization algorithm is presented which consists of cyclically heating and quenching by 
Metropolis and local search procedures, respectively. It works particularly well when it is applied to 
an archive of samples instead of to a single one. We demonstrate for the traveling salesman problem 
that this algorithm is far more efficient than usual simulated annealing; our implementation can 
compete concerning speed with recent, very fast genetic local search algorithms. 
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Combinatorial optimization problems have been a 
challenge for a long time. One of the most popular is the 
traveling salesman problem (TSP): how to find the short- 
est roundtrip through a given set of cities. It serves as 
a standard test for new optimization algorithms, for sur- 
veys see 00 Many combinatorial optimization tasks are 
of considerable practical importance. Thus algorithms 
are highly desirable which yield good approximations of 
the exact solution within a reasonable CPU time, and 
which need only a modest effort in programming. As 
such, simulated annealing [3J, |^ has become a standard 
procedure. Genetic algorithms have been increasingly 
successful, in particular such algorithms where local min- 
ima are the species considered, see 0,0,0 and references 
therein. Furthermore, e.g. alternately optimizing the sys- 
tem itself and various subsystems of it proved to be effec- 
tive 0| . Still, the danger of getting trapped in any local 
instead of the global minimum increases with the size of 
the problem and cannot be overcome by these heuristic 
methods. 

Simultaneously, exact methods have been developed 
further mainly using branch- and-bound and branch-and- 
cut ideas |3, llOj . Even a specific TSP instance includ- 
ing 7397 cities could be solved However, for a 
given size, the effort necessary to find the exact solu- 
tion can vary enormously from problem to problem; the 
TSPLIB95 11] includes an instance of 1577 cities which 
has not been solved exactly yet. The ideas basic in ex- 
act algorithms can also be used as background for fast 
local search procedures. Applied to randomly chosen 
states, such restricted branch-and-bound search imple- 
mentations can compete with simulated annealing with 
optimized schedule |12j . 

Here, we present an algorithm which combines the ad- 
vantages of simulated annealing with those of fast local 
search procedures. Additionally, our algorithm utilizes 
an archive of local minima similarly to a recently pub- 
lished hybrid-Metropolis procedure [l^. First, we ex- 
plain the general ideas of this algorithm. Then we illus- 
trate its efficiency by considering several TSP instances. 

Simulated annealing 0] can be understood as a ran- 



dom journey of the sample (i.e. the approximate solu- 
tion) through a hilly landscape formed by the states of 
its configuration space. The altitude, in the sense of a 
potential energy, corresponds to the quantity to be opti- 
mized. In the course of the journey, the altitude region 
accessible with a certain probability within a given num- 
ber of steps shrinks gradually due to the decrease of the 
temperature in the Metropolis simulation |14| involved. 
The accessible area, i.e., the correponding configuration 
space volume, thus shrinks until the sample gets trapped 
in one of the local minima. 

Deep valleys attract the sample mainly by their area. 
However, it is tempting to make use of their depth. For 
that, we substitute the slow cooling down by a cyclic pro- 
cess: First, starting from the lowest state obtained so far, 
we randomly deposit energy into the sample by means of 
a Metropolis process with a certain temperature T, which 
is terminated, however, after a small number of steps. 
This part is referred to as heating. Then we quench the 
sample by means of a local search algorithm. This is an 
iterative process which consists in searching systemati- 
cally the neighborhood (defined via the move class) of 
the actual state for states of lower energy, where always 
the first such state found is accepted. It stops when no 
lower neighboring state exists. Heating and quenching 
are cyclically repeated; the amount of energy (controlled 
by T) deposited in a cycle decreases gradually. This pro- 
cess continues until, within a 'reasonable' CPU time, no 
further improvement can be found. 

Our procedure has some resemblance to genetic local 
search algorithms I Mp l Bl as well as to the iterated Lin- 
Kernighan concept |li ll5l . But there are substantial 
differences. In genetic local search algorithms, crossover 
is the process to produce long jumps overcoming barriers 
in the configuration space. In the iterated Lin-Kernighan 
approach, the sample is kicked out of the local minimum 
by moves chosen from a specific class not considered in 
local optimization. Our algorithm, however, is not based 
on introducing a class of particular elementary moves for 
crossing barriers. Instead, it makes use of an inherent 
feature of combinatorial optimization problems: whether 
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or not a certain move is permitted, or reduces energy, de- 
pends on the actual state. Loosely speaking, the moves 
do not commute. Thus, provided T is large enough, even 
exciting the sample by a sequence of several rather prim- 
itive moves of a type used also in quenching should be 
sufficient for the cycle returning to its initial state only 
very seldom. This simple concept should be applicable 
to all problems which can be treated by simulated an- 
nealing. Moreover, this approach permits to make good 
use of certain complex moves, which would diminish the 
efficiency of heating due to their huge number, by con- 
sidering them only in quenching. 

When our algorithm is applied to a specific problem, 
several choices have to be made. The first one concerns 
the starting temperature, cf. [lTj. One can relate it to the 
energy difference between random and quenched states, 
or to the local misfit in quenched states, in the sense of 
frustration of a spin glass. 

The next issue is the amount of energy deposited into 
the sample while heating. Two contradicting demands 
have to be met: the gains of the previous cycles should 
be retained, but the modifications must be sufficiently 
large, so that another valley can be reached. Thus the 
heating process has to be terminated in an early stage of 
the equilibration. An effective method is to stop it after 
a fixed number of successful Metropolis steps. 

Our recipe for controlling T makes use of the rate of 
finding final states of the heating-quenching cycles which 
have lower energies than all previous ones: we keep T 
constant as long as this rate exceeds a certain value. The 
lower this value, the better should on average be the final 
approximate solution. 

The number of modifications performed within one 
heating process is small; moreover, only few of the cy- 
cles yield an improvement of the best state so far. Thus 
a considerable accclaration can be achieved by a reduc- 
tion of the move class considered in heating: after each 
change of T, we tabulate those moves which excite the 
best state so far by energies not larger than T. While 
heating, we consider only these moves. However, during 
quenching, we consider the related complete move class. 

Finally, one has to determine a stopping criterion. The 
most natural way is to compare the probability that the 
cycle returns the sample to its initial state (or to one with 
the same energy) with the probability of finding states of 
lower energy. 

We improve now the above approach by searching the 
configuration space starting the cycles at random from 
one of JV a local minima held in an archive , instead of 
from the best state so far. In this simultaneous search, 
the above described table of moves for heating comprises 
the exciations of all archive states by energies up to T. 
We reduce this table utilizing 'collective experience': If 
all archive states have a certain common feature, the 
global minimum is likely to also have this feature. There- 
fore, after changing T, we determine those degrees of free- 



dom which are frozen out, and eliminate all excitations 
from the table which concern one of those. This reduc- 
tion corresponds to the use of 'don't look bits' pj], or to 
the search for backbones 0] . 

In such a simultaneous search, it is tempting to sim- 
ulate evolution by mutation and competition: If the en- 
ergy of the final state of a cycle is smaller than that of 
the highest archive state, one would substitute this final 
state for the highest archive state. However, we observed 
that it is mostly preferable to replace the initial state of 
the actual cycle instead of the highest archive state. 

Concluding the general part we mention two other in- 
terpretations of such an algorithm. In the biological con- 
text, our procedure is of the genetic local search type 
\ll ki LZj - But it is based on mutations only, crossovers 
are not used. From an economist's viewpoint ft 
increases the performance by means of risk reduction on 
two levels: (i) After leaving a local minimum, we perform 
only a small number of modifications before we decide 
whether or not the continuation of this way is promising, 
(ii) Several valleys are investigated simultaneously. 

We have applied the above ideas to a series of sym- 
metric traveling salesman problems, as well as to the 
Coulomb glass. In both cases, the corresponding codes 
work very efficiently. Here, we focus onto the TSP, where 
the length of the roundtrip has the meaning of the en- 
ergy to be minimized. We use the following notions: tour 
and subtour denote closed roundtrips through all cities 
or part of the cities, respectively, whereas chains and sub- 
chains stand for tours and subtours with one connection 
eliminated, respectively. 

The efficiency of our procedure, as well as that of sim- 
ulated annealing, depends to a large extent on the move 
class taken into account. We consider two elementary 
heating processes: (i) cutting out from the tour a sub- 
chain, reversing its direction, and inserting it again at 
the same place; (ii) shifting a city from one to another 
position in the tour. However, we do not consider all 
such moves, but use the above introduced restrictions. 

In quenching, we choose between four possibilities con- 
cerning the kind of metastability to be reached: 

(a) stable with respect to reverse of a subchain, as well as 
to a shift of a city (unrestricted move class of heating); 

(b) same as (a), and stable with respect to cutting three 
connections of the tour, and concatenating the three sub- 
chains in a new manner; 

(c) same as (b), and stable with respect to first cutting 
the tour twice and forming two separated subtours, and 
connecting then these subtours after cutting two other 
connections; 

(d) same as (c), and stable concerning a restricted Lin- 
Kernighan search |2jj which consists in cutting the tour 
once, then several times alternately cutting the chain and 
concatenating the subchains, and finally connecting the 
ends of the chain again, where the number of trials to 
modify the chain is restricted to 1000. 
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The efficiency of our local search program rests on trying 
new connections according to increasing length, where 
appropriate bounds are utilized to terminate the search 
as soon as it becomes useless. Moreover, for stage (c), 
we demand that the sum of the lengths of the two closed 
subtours is shorter than the length of the original tour. 

Our algorithm includes several adjustable parameters. 
They all influence simultaneously both the quality of the 
final result and the CPU time needed. Nevertheless, we 
have observed that they affect only rather weakly the 
performance of the algorithm, i.e., the quality of the av- 
eraged final result obtained for the archive size iV a ad- 
justed so that the program finishes within a certain CPU 
time. For the TSP instances considered, the following 
values were good choices. First, we quench 50iV a times 
a randomly chosen tour, and set up the initial archive 
from the shortest N a of the tours obtained. As start- 
ing temperature, we choose the reduction of the length 
when quenching a random tour, divided by the number 
of cities. In heating, we modify the tour always 50 times. 
At each T, we perform 5N a heating-quenching cycles. If, 
during these cycles, an archive state has been replaced, 
we perform again 5N a cycles, and so on. Otherwise, T 
is reduced by a factor of 0.9. We terminate the process 
when it has returned 10iV a times to an archive state (or 
to a state with equal energy) without any substitution of 
an archive state. 

We tested our program considering five symmetric TSP 
instances out of the TSPLIB95 11], including between 
442 and 3795 cities. We performed several tests with 
different numbers of archive states, quenching always to 
metastablity (d) defined above. The results are given in 
Table I. They demonstrate the high efficiency of our al- 
gorithm: For the instances with 442, 532, and 783 cities, 
the known exact optima were reproduced. The instances 
with 1577 and 3795 cities have not been solved exactly 
yet. Here, our best tour lengths agree with the known 
upper bounds for the optimum tour length [llj . 

However, our algorithm yields generally only an ap- 
proximate solution. Thus the truly important property 
is the relation between mean quality of the solution and 
the necessary CPU time. We compared the performance 
of our program with that of the genetic local search al- 
gorithm |jj by Merz and Freisleben - a significantly im- 
proved version of their winning algorithm of the First In- 
ternational Contest on Evolutionary Optimization |2l| -, 
and with the performance of the iterated Lin-Kernighan 
code by Johnson and McGeoch, illustrated by Table 16 
ofQ (Our CPU is by factors of roughly 1.7 and 3 faster 
than the CPUs used in0 andfl respectively.) 
- Compared toQ, our code has roughly the same efficiency 
for the three instances with several hundred cities. How- 
ever, for the instances with 1577 and 3795 cities, our pro- 
gram is more efficient. In particular, for the task A3795, 
it yields a significantly better mean value in less than one 
fourth of the CPU time. 



Problem iV a 


J^min 


^max 


^mcan 


TCPU 


pcb442 


1 


50778 


51004 


50849 


25 


pcb442 


3 


50778 


50912 


50794 


69 


pcb442 


5 


50778 


50778 


50778 


145 


att532 


1 


27686 


27778 


27722 


32 


att532 


3 


27686 


27732 


27707 


88 


att532 


5 


27686 


27715 


27694 


191 


att532 


8 


27686 


27705 


27692 


372 


att532 


12 


27686 


27693 


27686.4 


645 


rat783 


1 


8809 


8829 


8818 


38 


rat783 


3 


8806 


8823 


8812 


101 


rat783 


5 


8806 


8810 


8806.7 


213 


rat783 


8 


8806 


8809 


8806.1 


427 


rat783 


12 


8806 


8806 


8806 


759 


A1577 


1 


22254 


22289 


22264 


251 


111577 


3 


22250 


22267 


22256 


808 


H1577 


5 


22249 


22261 


22252.5 


1730 


H1577 


8 


22249 


22255 


22251.8 


3350 


H1577 


12 


22249 


22253 


22250.4 


6070 


H3795 


1 


28775 


28890 


28809 


2100 


A3795 


3 


28775 


28830 


28790 


6300 


H3795 


5 


28772 


28784 


28776.5 


13500 


H3795 


8 


28772 


28780 


28774.8 


25900 



TABLE I: Dependence of the tour length of the approximate 
solution, and of the CPU time on the archive size iV a for 
five instances out of the TSPLIB95 [Till . Smallest, largest, 
and mean tour lengths, L m i n , L max , and L moan , are given for 
series of 20 runs, tcpu denotes the CPU time in seconds for 
a single run using one PA8000 180 MHz processor of an HP 
K460. The local search guaranteed stability with respect to 
demand (d) defined in the text. 

- Compared to Q, our program is by a factor 5 to 10 
slower for pcb442 and att532 if the accuracy of our re- 
sults for N a = 3 is achieved. The performance gap seems 
to shrink with increasing accuracy demand |22|. How- 
ever, for A3795, our code performs considerably better. 
The clear advantage of our program for A3795 in both 
comparisons can arise from particular robustness, or 
from good scaling properties. Further investigations are 
needed to clarify this point, and to decide whether the 
advantage is caused by the use of thermal cycling to over- 
come barriers, or by a feature of our local search proce- 
dure. 

One should also make comparison with the state of 
the art exact algorithms. For the Padberg-Rinaldi 532 
city problem, the branch-and-cut program by Thienel 
and Naddef, one of the presently fastest exact codes, 
needs 16.5 minutes at a SPARC10 '23\, corresponding 
to roughly 4 minutes for our CPU. Utilizing an archive 
of 12 states and cyclically quenching to stage (d) defined 
above, we performed 100 runs. Our Monte Carlo ap- 
proach reproduced the optimum tour length 27686 in 95 
of these runs, where it needed 11 CPU minutes for one of 
them. In the residual five cases, we obtained tours with 
length 27693. 

We have shown above that thermal cycling can be used 
as basis for state of the art Monte Carlo optimizations. 
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FIG. 1: Relation between CPU time, tcpu (in seconds for 
one PA8000 180 MHz processor of an HP K460), and devi- 
ation, 8L — Z/ m0 an — 27686, of the obtained mean approxi- 
mate solution from the optimum tour length for the Padberg- 
Rinaldi 532 city problem. □: repeated quench to stability 
with respect to demand (a) defined in the text; A: simulated 
annealing; x, Vi +> * : thermal cycling with archives of vari- 
ous size, and local search concerning conditions (a), (b), (c), 
and (d), respectively. In all cases, averages were taken from 
20 runs. Errors (lcr-region) are presented if they exceed the 
symbol size. The straight lines are guides to the eye only. 



Now, we demonstrate that this physical approach can 
be far more effective than standard concepts, even if 
only a comparatively small amount of time is invested 
in programming. For the Padberg-Rinaldi 532 city prob- 
lem, Fig. 1 confronts four versions of our algorithm, dif- 
fering from each other by the complexity of the local 
search, with usual simulated annealing, and with repeat- 
edly quenching randomly chosen states. (Simulated an- 
nealing uses the same moves as heating above, but trials 
were restricted to the 30 nearest cities.) The figure shows 
that the thermal cycling procedure proposed is far faster 
than usual simulated annealing, even if the local search 
leads only to the 'primitive' metastablity (a). 

Concluding, we have shown that simulation of thermal 
cycling is an effective tool for treating combinatorial op- 
timization problems. Its efficiency in finding very good 
approximations of the optimal solution was illustrated by 
considering the traveling salesman problem. Our algo- 
rithm, which is well appropriate for parallelising, should 
also be useful for a broad class of other combinatorial 
optimization problems. 

This work was supported by the SMWK and DFG 
(SFB 393). We are particularly indebted to B. Freisleben 
and S. Thienel for valuable literature hints, as well as to 
M. Pollak and U. Rofiler for a series of critical remarks. 
Moreover, discussions with F.-M. Dittes, M. Golden, 
D.S. Johnson, S. Kobe, A. Hartwig, M. Ortuno, 
M. Richter, and J. Talamantes were very useful. 
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